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Abstract 

The structural flexibility of low dimensional nanomaterials offers unique opportunities for study¬ 
ing the impact of strain on their physical properties and for developing innovative devices utilizing 
strain engineering. A key towards such goals is a device platform which allows the independent 
tuning and reliable calibration of the strain. Here we report the fabrication and characterization of 
graphene nanoelectromechanical resonators(GNEMRs) on flexible substrates. Combining substrate 
bending and electrostatic gating, we achieve the independent tuning of the strain and sagging in 
graphene and explore the nonlinear dynamics over a wide parameter space. Analytical and nu¬ 
merical studies of a continuum mechanics model, including the competing higher order nonlinear 
terms, reveal a comprehensive nonlinear dynamics phase diagram, which quantitatively explains 
the complex behaviors of GNEMRs. 
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Understanding and controlling the effect of strain on the physical properties of low dimen¬ 
sional nanomaterials has been a topic of great interest in recent years. Various types of strain 
have been theoretically predicted to exert fundamental impact on the electronic properties, 


including flexural phonon scattering[l|, gauge field formation|2| and modification of lattice 
symmetry Through strain engineering, the mechanical, electronic, and optical properties 


of a nanomaterial may be tailored to better suit its functionality. While extensive work has 
been carried out to study those properties via optical QB and local probe techniques!^ j^, 
it is ultimately desirable to achieve strain engineering in electronic devices}^. The major 


hurdles in achieving this goal are the absence of an independent tuning knob for strain in 
an electronic device and a reliable approach for characterizing the strain. 

A promising device platform for realizing strain engineering is a suspended field effect 
transistor (FET)-like nanoelectromechanical resonator(NEMR). Graphene, for example, has 


been intensively s 


elastic modulus 


ndied in NEMR devicesjo] 1^, its large mechanical strength ll|, high 
12j and low mass density. With graphene suspended from the substrate. 


these devices potentially allow probing its intrinsic properties. In addition, the strain in 
graphene can be tuned electrostatically through the gate bias. And the magnitude of the 
strain can be estimated from the gate dependence of the resonant frequency of graphene. . 

Despite these advantages, a major limitation of tuning strain through electrostatic force 
is that it unavoidably affects charge doping and hence the electronic properties of graphene. 
Moreover, the electrostatic strain is always tensile. The range of such strain tuning is also 
hmited(~0.01% at 10v{l^), because large gate bias can cause collapse of graphene channel. 
This tuning range is usually much smaller than the residual strain from the fabrication 
process (~0.1%jl5|). 

Electrostatic tuning can also complicate the analysis of strain by inducing stron g no nlin- 
earity under large gate bias. Typically, large deformation 1^ , clamping geometry 17| and 


external force 


18| introduce higher order terms of displacement and restoring force. These 


terms play significant roles in the equation of motion and give rise to deviation from simple 
harmonic behavior. In FET-like GNEMRs, tuning the strain by perpendicular electrostatic 
force inevitably causes the sagging of graphene, which affects the resonance dynamics. Un¬ 
der large driving amplitude, broadened nonlinear resonance has been observed. Depending 
on the random residual strain and sagging, the broadened resonance shifts either to higher 


frequencies (i.e. hardening behavior) 


13| or to lower frequencies (softening behavior) 191 
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FIG. 1. Device fabrication and strain tuning, (a) Process flow diagram for device fabrication, 
(b) False color scanning electron micrope (SEM) image of a suspended graphene sheet with Au 
contacts and bottom gate, (c) Mechanism for strain tuning by bending the substrate, h is the 
deflection of substrate, which is negative for compressive surface strain, positive for tensile strain. 


with increasing oscillation amplitude. A comprehensive study of the relation between the 
strain and nonlinear dynamics is therefore important for determining strain in the nonlinear 
regime. However, further experimental study has been restricted, due to the lack of an 
effective tuning knob for the strain j^. 

Here we develop a device platform for strain engineering, by fabricating FET-like local¬ 
gated GNEMRs on flexible substrates. Bending the substrate enables the independent tuning 
of the strain. By analyzing the gate bias and substrate bending dependence of the nonlinear 
dynamics, the strain in graphene is reliably determined. 

The flexible substrate used here is 125 fim thick Polyimide, which maintains flexibility 
over a wide temperature range (~4K - 650K). The fabrication of a GNEMR on the Polyimide 
substrate is depicted in Fig. HJa). First, we evaporate 6 nm of Ti onto the Polyimide 
substrate to avoid the static charging in electron beam lithography (EBL). After patterning 
the substrate with 1pm wide, 30 nm thick gold gate electrodes by EBL, we coat the substrate 
with 300 nm polymer spacer (PMMA). Simultaneously, a polydimethylsiloxane/transparent 
tape/methylmethacrylate(MMA) stack is prepared on a transparent glass slide 2l|, with 
graphene flakes exfoliated onto MMA. Using a micro-manipulator, the MMA membrane 
with graphene flakes is transferred onto the PMMA covered substrate, with the chosen 
graphene flake aligned to the gate electrode. Then following an etching-free suspending 
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procedure 22j, we fabricate the electrical contact pads using EBL. Part of the polymer layers 
(both MMA and PMMA) are cross-linked so that the graphene is strongly clamped at the 


contacts 


23] 241. Following metalization with Cr(l nm)/Au(200 nm), the standard “lift-off’ 


process is carried out using acetone, after which the graphene is supported by the 3D metal 
contact “pillboxes”. In the end the exposed thin Ti conducting layer is etched away using 
diluted HP acid. Throughout the whole lift-off and etching process, the device goes through 
the following solvents successively: acetone, isopropyl alcohol(IPA), de-ionized(DI) water, 
diluted HP, DI water and IPA. Finally, it’s taken directly out of hot Hexane. From the 
thickness of the polymer spacer and the gate electrode, the distance between graphene and 
the gate d is estimated to be 270 nm. Fig. [I](b) shows the SEM image of a hnished device 
where graphene is suspended over the substrate by the 3D structure of the metal contacts, 
clamped by cross-linked polymer. In all our devices, we select rectangular graphene flakes 
with sharp and straight edges to avoid edge modes and local buckhng[25|. Bending the 
substrate with a deflection h (positive for outward bending and negative for inward bending), 
the surface strain of the substrate modihes the source-drain distance, hence tuning strain in 
graphene (Fig. [T](c)). 

For resonance measurements, we have adopted the electrical mixing method reported in 


a previous work 


13| . In this method, the oscillation is characterized by the mixed current 


through graphene (excluding background), which is proportional to the oscillation amplitude 
of the graphene membrane. The measurement is carried out in vacuum below 50 mTorr at 
room temperature. The deflection of the substrate h is controlled by adjusting a plunger 
with a threaded rod, which is coupled to the ambient through a vacuum rotary feed-through. 


The detailed setup and the circuit diagram are shown in the Supplementary Material 


26 |. 


All measurement results presented in this paper are obtained from a device with graphene 
channel of width W = 2 /rm and length L = 1 /rm. 

We hrst study the resonance of the GNEMR device without bending the substrate. 
Fig.EKa) shows the mixed current as a function of frequency with gate bias Vq tuned from 
-4 V to -10 V. The corresponding resonant frequency fr increases monotonically from 108 
MHz to 122 MHz. Sensitivity of fr on Vq indicates that the observed resonance is from the 
graphene and not from the electrodes (which is almost Vq - independent). Both the large 
values of fr and the lineshape of hardening behavior imply signihcant tensile residual strain 
in graphene. The resonant frequencies for various gate bias values are shown in Fig. [2](b) 
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FIG. 2. Dependence of resonant freqnency on gate bias, (a) Mixed current (proportional 
to the oscillation amplitude) versus frequency at various DC gate bias Vq. As the bias voltage 
changes from -4V (black) to -lOV (yellow), the resonant frequency increases monotonically. (b) 
Experiment data and simulation results of resonant frequency versus gate bias at various values of 

A. 


(see the data labeled by squares). The values are symmetric around Vq = 0 and qualita¬ 
tively consistent with previous report [^. For the device studied here, since the hole branch 
shows higher transconductance and gives better resonance signal than the electron branch 
(see Fig. S 1(b) in the Supplementary Material 2q), we focus on the measurement at negative 
gate biases. 

Next we focus on how the resonance depends on the strain A. Here we define A as the 
strain in graphene at zero gate bias, which takes into account both the residual tension 
and the substrate bending. The tensile residual strain indicates A > 0 without bending. 
Fig.EKa) shows the mixed current as a function of the frequency at Vq = —7V, as the tensile 
strain is reduced by inward bending of the substrate. Here the resonant frequency decreases 
from 114 MHz to 86 MHz as |h| increases to 159 /um. Further inward bending results in the 
increasing of fr up to 92 MHz. This non-monotonic dependence of fr on A generally holds 
for different Vq, as summarized in Fig. 111(b). Moreover, as Vq is lowered, A corresponding to 
the minimum resonant frequency gets smaller. Here A is quantitatively calibrated through 
numerical simulation discussed later. We note that the observed A and Vq dependence of 
the resonant frequency is reproducible upon repeated inward bending and releasing of the 
substrate, indicating a good mechanical clamping of the graphene channel and an absence 
of structural slippage. 

In conjunction with the non-monotonic behavior of the resonant frequency, the resonance 
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FIG. 3. Dependence of resonant frequency and lineshape on strain A. Different color 
corresponds to different A. (a) Experiment results. The direction of the arrows indicates increasing 
inward bending (i.e larger \h\ and smaller A), (b) Simulation results. The simulated frequency and 
lineshape of resonance peaks match with the experiment results. Both experiment and simulation 
results show non-monotonic dependency of resonant frequency and transition from hardening to 
softening resonance when A is decreased close to zero. 


lineshape also undergoes a transition around the minimum resonant frequencies in Fig. [3](a). 
Without substrate bending, the resonance shows a hardening behavior. As \h\ increases 
to 136 fim, the resonance peak becomes symmetric. At the minimum fr, the softening 
resonance emerges. The softening behavior gets increasingly pronounced with further inward 
bending, as indicated by the steeper rise of the resonance amplitude at the tippin g p oint. 


While softening and hardening behaviors have been randomly observed in GNEMRs[l^ 
here we achieve hne control of the nonlinear dynamics by combining lateral stretching and 
vertical electrostatic pulling. 

To quantitatively understand the observed resonance behavior, we model graphene as a 


continuum membrane with ID approximation 


27] 28] 29] 30], which is illustrated in Fig.[5](a) 


We consider the driving gate voltage to be V (t) = Vo + 6V cos {cot), under which the deflection 
of graphene is approximately a combination of a parabolic static deflection and an oscillating 
fundamental mode: z{x,t) = zo{l — /L'^) + u{t)cos{7ix/L). Here zo is the displacement 

at the center of graphene, L is the source-drain distance. Assuming the strain-free length of 
the graphene beam is Iq, the strain A = (L — lo)/lo- We then obtain the following nonlinear 
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FIG. 4. Dependence of resonant frequency on gate bias and A. Solid lines are fits from 
model. Symbols are experiment data, (a) Linear fitting of versus . (b) Resonant 

frequency fr versus strain A at various gate biases. 


differential equation for the oscillatory component u: 

du , . 9 ^ AenVodVcosiujt) 

PysTT + r]^ + ku + Au^ + Bu^ + ^ = 0 . 

oB ot Tid^ 


( 1 ) 


Here p is the 2D mass density, p characterizes the damping strength of the resonator, (i(s> z) 
is the vertical distance between the source-drain and the gate electrode, k = + 

(8/3-|-512/7r"‘)^QL“^], A = 167rii^2;o-h~'^ and B = /4)EL~^, in which E is the 2D Young’s 

modulus of graphene. Note that these expressions assume, as appropriate for the geometry 
of our device, that zq and u in the graphene deflection are much smaller than the distance d 
to the gate electrode. In this regime, electrostatic energy of the graphene-gate capacitance 
gives only negligible contribution to the coefficients fc. A, B. In particular, for the effective 
spring constant k, the electrostatic contribution: eoVod~^ = 16zo/dEL~^[X + 8/3 zqL~‘^], is 
of the order oi zo/d I'm. relative to the elastic contribution. If ~ d, electrostatic 
energy can change the behavior of the resonance Ml] ■ 

The main resonant frequency can be characterized by fr = 1/2 tt^J~ kjp. To qualitatively 
understand how fr varies with the strain A, we can consider the dependence of the effective 
stiffness of the graphene channel on the deformation along the z-axis, described by an ef¬ 
fective spring constant fee// — d{Td‘^z/dx‘^)/dz)\x=o- The presence of a minimum frequency 
can be understood considering the competition between the tension T and the curvature 
d‘^z/dx^. When A > 0 and is large, the sagging is small and the oscillatory motion is per¬ 
pendicular to the direction of tension. The behavior of the resonator is similar to that of a 
classical string whose resonant frequency decreases with decreasing tension. On the other 
hand, when A gets smaller, the sagging gets larger, increasing the curvature of graphene. 
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FIG. 5. Continuum mechanics model for doubly clamped graphene resonator, (a) 

Schematic of the graphene resonator, zq is the static displacement of the center of graphene. L 
is the distance between S-D electrode clamps. The shape of graphene is described by z(x). d is 
the distance between S-D and gate. (b) Nonlinear contribution to the potential energy plotted as 
a function of average oscillatory displacement u. (c) Nonlinear dynamics phase diagram. Three 
characteristic regimes - hardening (H), softening(S) and the non-monotonic softening (WS) regime 
- are separated by the white solid lines on the diagram. The black dashed line highlights the 
frequency minima at fixed gate biases and is close to the hardening-to-softening transition line. 
The grey dashed lines correspond to the frequency minima at fixed A. 


This significantly conples the oscillatory motion with the longitndinal deformation. As a 
resnlt, fr starts to increase with decreasing A. A more rigorons analysis taking into ac- 
connt the oscillatory component (see the Snpplementary Material 2^ ) hnds the minimnm 
freqnency and the corresponding A to be: 


_ 1 4 1/6 

jrain r>_ 1 r> r 4 ..74 ^0 / > 


27r ^ 2p3L4d4 

\ — 2.1/3, 141^x2/3 


( 2 ) 


r3 9 ^ rT/2 

^TT 1922 ^ 2Ed‘^ 

The transition from hardening to softening behavior with decreasing A can be nnder- 
stood qnalitatively considering the competition between the qnadratic and the cnbic term 
in Eq.([T]). Under large tension, zg — 0 and the qnadratic term is negligible. So the ma¬ 
jor contribntion to the “mean” k^ff is from the positive cnbic term, which increases with 
increasing amplitnde. This gives rise to the hardening Dnffing oscillator behavior. On the 
other hand, nnder small tension, the presence of a gate bias canses significant sagging, hence 
the qnadratic term dominates. This term induces a tilt in the potential energy, as illustrated 
in Fig. |5](b). As a result, with increasing oscillatory amplitude, the average oscillatory dis¬ 
placement u shifts towards more positive values, where the curvature of the potential well 
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(the keff) is smaller. This gives rise to a decreasing resonant frequency with increasing 


be 


lavior. Following a more detailed calculation shown in the 


amplitude, hence softening 

Supplementary Material 2^ j32| , the switching between the two behaviors, in the limit of 
small driving force, happens at A = This shows similar depen¬ 

dence on Vo as Eq.(l2]), but with a slightly larger prefactor. 

To obtain a quantitative comparison between the model and the data, we hrst evaluate 
the mass density of the sample. In agreement with Eq.(l2]), the data of versus 
follows a linear dependence which extrapolates to the origin(Fig. 111(a)). The slope, which 
is proportional to 1/p, provides an estimation to the mass density of the channel. Adopt¬ 
ing E^d = 1 TPa llj] and using the device geometry d = 270nm and L = l/rm, we obtain 
the mass density to be 5po, in which po = 7-4 x is the mass density of mono- 

layer CTa^ene. The large mass density, which has been frequently observed in previous 


work 131] 19|, can be attributed to the resist residue and other contaminants. Using the 


estimated mass and Eq.([T]), we are able to numerically simulate the gate bias dependence 
of the resonant frequency at various A. For the simulations, we use a dissipation factor of 
p = 50 and the experimental driving amplitude 5V = 300mU, but vary the value of A as 
a htting parameter. The simulation results agree well with the data from experiment, as 
shown in Fig. [2](b). We also compare the simulated line shapes of the resonance peaks with 
the measurement results. As shown in Fig. l3](b), the model reproduces our observation in 
Fig. EKa) remarkably well. The suspended graphene device analyzed here shows a residual 
strain of ~0.042%, which is relaxed down to ~0% through substrate bending. The small 
discrepancy between the simulation and the experiment results may be attributed to the 
errors in the device geometry estimation and the simplihcation of our model, which neglects 
the nonlinear damping [l^ and the higher order modes. 

A comprehensive parametric dependence of the nonlinear dynamics in the GNEMRs is 
shown in Fig. [5](c). Based on numerical calculations, the color-coded resonant frequency 
is plotted as a function of gate bias and A. Three characteristic regimes can be identihed 
on the diagram. As observed in our devices, the hardening(H)-to-softening(S) transition 
happens close to the minimum frequency line (dashed black line). The regime, not covered 
in this work but has been previously reported 19|, is the WS regime for negative A. Here 
the gate-dependence of resonant frequency become non-monotonic with a “W” shape. The 
nonlinear dynamics phase diagram covers the so-far observed behaviors of GNEMRs and 
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may serve as a useful guide to the analysis of ID-like mechanical resonators in general. 

In summary, we have demonstrated a FET-like suspended graphene nanoelectromechani¬ 
cal resonator on a flexible substrate, where the strain can be independently tuned by bending 
the substrate. This design allows for the observation of nonlinear dynamics in GNEMRs, 
including non-monotonic strain-dependent resonant frequency and the transition from hard¬ 
ening to softening resonance behavior. The device platform developed here allows the in¬ 
dependent tuning of strain for measurements on the electronic properties of the channel 
materials, without electric held doping effects. This vastly broadens the scope of the strain¬ 
engineering of 2D nanomaterials. 
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